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ABSTRACT 

We construct a family of semi-analytic models for young open clusters, including 
a gaseous component and varying assumptions about the distribution function for the 
stellar component. The parameters of these models are informed by observed open 
clusters and general theoretical considerations regarding cluster formation. We use this 
framework to estimate the fraction J^i, of the stellar component that remains gravita- 
tionally bound after the gaseous component disperses. The remaining stellar fraction J^^ 
is a smooth function of the star formation efficiency e*, and depends on the distribution 
function of the stars. We calculate the function .^^(e,,,) for this class of open cluster 
models and provide fitting formulae for representative cases. 

Subject headings: open clusters and associations: general - stars: formation 

1. INTRODUCTION 

In this paper, we study the equilibrium structure and early evolution of young open star 
clusters. Our long term goal is to develop a unified treatment that includes cluster formation, 
removal of the gaseous component, and the longer term evolution of the system. A unified picture 
of young cluster evolution can be useful in several different contexts. On one hand, we can consider 
an open cluster as an astrophysical entity and study its birth, evolution, and ultimate demise. On 
the other hand, we can study the effects of the cluster setting on the star formation process. In this 
present work, we construct models for young clusters in which a gaseous component is still present. 
We then study how the clusters react to gas removal. Our results are applicable to robust clusters 
that live for relatively long times (~ 100 Myr); perhaps 10% of all star formation takes place in 
such robust cluster environments (see our companion paper Adams &: Myers 2000). 

A collection of previous papers have addressed the question of whether or not a cluster can 
remain gravitationally bound after the removal of its gaseous component (e.g.. Hills 1980; Mathieu 
1983; Elmegreen 1983; Lada, Margulis, &: Dearborn 1983). Most previous work uses an approach 
based on the virial theorem, however, and does not explicitly include the distribution function for 
the stars. In the case of rapid gas removal and a cluster that begins in a state of exact virial 
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equilibrium, the cluster expands by a factor /ex after the gas is removed. In the simplest models, 
this factor /ex is given by 

where is the mass fraction of the stellar component in the original cluster. In this approximation, 
the cluster remains bound for e* > 1/2; as e* 1/2, /cx oo and the cluster becomes unbound. 
This approach has been generalized to include varying assumptions about the speed of gas removal 
and departures from exact virial equilibrium, e.g., stellar speeds that are less then the virial speeds 
in the cluster potential (see, e.g., Verschueren 1990; Verschueren &; David 1989; Elmegreen & 
Clemens 1985; Lada et al. 1983). In this present work, we take into account the distribution 
function for the stars and, separately, the density distribution of the gaseous component. These 
complications lead to a much wider range of possible behavior. Because the cluster stars have a 
velocity distribution, the low speed stars in the tail of the distribution survive as a gravitationally 
bound entity even if < 1/2; the high velocity stars in the opposite tail escape even if > 1/2. 
For a given distribution function and a given star formation efficiency e*, a fraction J^*(e*) of the 
stars will thus remain bound after gas removal. The function varies smoothly with star 

formation efficiency rather than exhibiting singular behavior. 

This paper is organized as follows. In §2, we construct equilibrium models of young open 
clusters including standard forms for the stellar distribution function, a gaseous component, and 
anisotropy parameters. We study the effects of gas removal in §3; in particular, we calculate the 
fraction of stars remaining after gas leaves the system. We compare our results with observed 
clusters in §4 and then conclude, in §5, with a summary and brief discussion of our results. In the 
Appendix, we also present a crude cluster formation theory, which informs the theoretical models 
in the main text. Starting with the existing theory for the collapse of molecular cloud cores, we 
scale up the solutions to describe the collapse flow that forms a star cluster. Using this formalism, 
we estimate the anisotropy of the velocity distribution for forming clusters. 



2. EQUILIBRIA OF STAR CLUSTERS CONTAINING GAS 

In this section, we construct a class of cluster models that incorporate both stars and gas. In 
this standard approach, the structure of clusters is determined by two differential equations. The 
first is the Poisson equation for the gravitational potential 

V2$ = 47rGptot, (2) 

where the total density ptot includes both stars and gas. The second is the collisionless Boltzmann 
equation, which takes the form 

where / = /(x, v,t) is the distribution function for the stars. 
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2.1. Models with isotropic velocity dispersion 

As a starting point, we construct equilibrium models with an isotropic velocity dispersion. 
Making straightforward extensions of the general treatment outlined in Binney &; Tremaine (1987), 
we include a gaseous halo in the formalism. In particular, we use lowered isothermal cluster models 
by assuming a stellar distribution function /(e) of the form 

/(e)=pi(27ra2)-3/2[e^/-^-l], (4) 

where the relative energy e is defined by £ = ^' — and where the relative potential is related 
to the true gravitational potential $ through ^ = — <1> + $o- The constant pi sets the density 
scale and the parameter a sets the scale for the velocity dispersion (notice that a / (i;^)^/^). By 
assuming a distribution function which is a function of the energy only, we automatically obtain a 
solution to the collisionless Boltzmann equation [3] for equilibrium configurations with df / dt = 0. 
We have also implicitly assumed an isotropic distribution in velocity space - we consider the more 
general case of anisotropic models in a subsequent section. 

Prom the distribution function /(e), we can determine the stellar density in terms of the 
potential. Using the resulting density in the Poisson cqiiation, in conjunction with the additional 
terms resulting from gas, we can find the potential as a function of radius in the cluster, and then 
find the stellar density as a function of radius. For convenience, we define new variables 

^ = ^S/a\ ^ = r/ro, and ro = [9a'' /47rGpof\ (5) 

where po is the central density of the stellar component of the cluster. The variable ro is thus the 
King radius and plays the role of an effective core radius for the stellar component (see Binney & 
Tremaine 1987). 

Using the newly defined variables, we can evaluate the density of stars p*{r) in terms of the 
reduced potential ipir), i.e.. 



Pi 



e^Erf(v/V^) - 2VVV^(1 + ^V')] , (6) 



where Erf (2;) is the error function (e.g., Abramowitz &; Stegun 1972). The central density can thus 
be written 

PO = pi [e^°Erf ( V^) - 2^/^{l + ^^0)] , (7) 

where V'o = V'(O) = *(0)/(t^ determines the depth of the gravitational potential well at the cluster 
center. We find it convenient to define a function 

5*(x) = e^ Erf (Vi) -2^^(1 + 2x73), (8) 
so that = pig*iip) = pofl'*(V')/fl'*(V'o)- 
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To specify the gas contribution, we use a density distribution much hke an isothermal sphere; 
in particular, we consider a gaseous halo around the cluster with a density profile of the form 

a? 

where a is the effective sound speed of the gas and where rog is the core radius for the gas. For 
simplicity, wc assume that the gas and the stars have the same core radius ro; in principle, we could 
include an additional parameter CKg* = rog/ro in the analysis. 

The new version of the Poisson equation (which determines the cluster structure) takes the 

form 

where the constant A determines the contribution of the gas relative to that of the stars and is 
defined by 

We thus have a two parameter family of cluster models. To specify each model, we need to choose 
the value of '0o to set the depth of the gravitational potential well and the value of A to set the gas 
contribution. 

As we integrate equation [10] outwards, both the potential if) and the stellar density p* decrease 
and eventually become zero at some outer truncation radius r^. [Notice that the stellar density 
approaches the form p*{ip) ~ (S/lS-y/vr) tp^^'^ in the limit 0.] At this outer boundary, the 

true gravitational potential $ is given by ^{vt) = —GM{rT)/rT, where M{rT) is the total mass 
enclosed, including both stars and gas. The value of the true potential at the origin is then 
$(0) = ^{rr) — *(0), where *(0) = V'oO'^- As the input parameter i/jq grows larger, the outer 
truncation radius rr grows accordingly; the mass M{rT) also increases, and so does the overall 
depth of the potential well as determined by |$(0)|. 

A typical density profile is shown in Figure la (for a model using ipQ = 10 and A = 1/2). 
For this case, the mass in stars and the mass in gas are roughly comparable (specifically, 44% 
stars and 56% gas). The dashed curve shows the stellar component (eq. [6]) and the dotted curve 
shows the gaseous component (eq. [9]). The total density is given by the solid curve. After the 
gas is removed, the remaining stellar component contains ~ 73% of the starting inventory, with its 
resulting density distribution depicted by the dot-dashed curve (the formulation of this gas removal 
calculation is presented §3). For this class of models, the density in the cluster center is dominated 
by stars and the outside is dominated by gas; the star formation efficiency thus increases toward 
the inside. This structure is roughly consistent with that expected from a collapse model of cluster 
formation as outlined in the Appendix. Even though the density is strongly peaked toward the 
center, most of the mass in the cluster - in both stars and gas - resides in the outer portion of the 
cluster. Figure lb shows a collection of density profiles with varying depths of the gravitational 
potential, determined by tpQ which lies in the range 6 < -00 < 14. Similarly, Figure Ic shows a 
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collection of density profiles for ipo = 10 and varying amounts of gas, determined by A which Hes 
in the range < A < 2. 

In this present construction, the stars hve in a potential well determined by both their own 
self-gravity and by the gravity of the accompanying gas. The stars have a distribution function 
which was chosen to be a function of an integral of motion (the energy). Through the Jeans 
Theorem, the distribution function is thus a valid solution to the collisionless Boltzman equation 
which determines equilibrium configurations of stellar dynamics (see Binney & Tremaine 1987). 
Although physically motivated, the gaseous halo distribution [9] was an arbitrary choice and will 
not in general be in perfect hydrostatic equilibrium with the potential obtained in solving equation 
[10]. To address this potential discrepancy, we assume that support for the gaseous halo is given 
by additional pressure terms produced by external agents such as magnetic fields. This approach is 
adopted in the present version of the calculation. In principle, we could also solve self-consistently 
for the density distribution of the gas through an iterative procedure. 



2.2. Models with anisotropic velocity dispersion 

In this section, we consider the more realistic case of anisotropic models. For an anisotropic 
velocity dispersion tensor, the distribution function / is a function of both the energy e and the 
angular momentum L. To obtain a natural extension of the formulation presented in the previous 
section, we use a distribution function of the form 

/(£, L) = pi(27ra2)-3/2 giL^rW) {e'/'^' - 1) , (12) 

where the function g and the anisotropy radius rj^ determine the degree of departure from isotropy. 
In the limit rA oo, g ^ 1 and we recover the isotropic models of the previous section. In 
our cluster formation scenario (see the Appendix), the anisotropy radius va is determined by the 
centrifugal radius Rc of the collapse flow so that rA ~ Rc (where Rc is evaluated at the end of 
the cluster formation epoch). 

With a distribution function of the form [12], the stellar density can be written 

p^ = p^^±J siu-qdr] j v"^ dv g{r\'^ sin^ r}/r\) (e'^-^'/s _ i) ^ (13) 

where we have scaled out the a dependence. Motivated by our model of cluster formation from a 
collapse flow, we use an anisotropy function g{L'^) that decreases as one power of r (and hence L) 
for large radii. We thus choose the particular form 

^ " (l + L7r2)l/2 = (l + ^2^2gij,2^/^^)l/2 = (l + «2^2sin2^)V2 ' ^^^^ 

where we have defined an anisotropy parameter a = r/rA- With this choice for g{L^), the angular 
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integral becomes 

' ' - — atan{av) . (15) 



r7\ 
JO 



lo [l+a2i;2sin2?7]^/2 av 
The stellar density takes the form 

= puf^"-^ vdvatan{av){e'^-'''/^ - 1) = p^/^!^ , (16) 

\ r Jo V TT r 

where we have defined the integral Ip in the final equality. The integral Ip can be partially evaluated 
to obtain 



Ip=^-atan{a^/2^)[l + ^ + ^]+ae^Kp, (17) 



where the remaining integral Kp is given by 



p-«V2 



Although the integral Kp cannot be evaluated in terms of elementary functions, we can obtain 
a good approximation using asymptotic analysis (Bleistein &; Handelsman 1986). In the limit of 
small a, we invoke Laplace's method and keep only the leading order term to find the following 
analytic expression 

Kp « y|(l + 2a2)-i/2 Erf[VV'(l + 2a2)] . (19) 

This approximation becomes exact in the limit of small a (small departures from isotropy). In the 
opposite limit of large a, 

Kp K atan{ay/2ijj)/a . (20) 

The stellar density contribution is thus given by the combination of equations [16 - 20]. Since we 
have made an approximation in evaluating the density (equations [19, 20]), the distribution function 
corresponding to this density profile is not exactly of the form given by equation [12], although it 
is close and it has the correct behavior in the limiting regimes. 

The form of the Poisson equation, which determines the potential, now takes the form 

^ d^y^ dc K{^o) 1+e'' ^ ' 

where the function /i*(x) is defined by 

\/2x 1 r 1 1 

h^{x) = e^Kp{x, a) + atan(a-\/2x) 1 + x + . (22) 



2a2 a 



The parameter a = rjrA = ^.To/r^, so we define an overall anisotropy parameter /3 = tq/va- Our 
scenario for cluster formation indicates that both the anisotropy radius ta and the cluster core 
radius ro are approximately given by the centrifugal radius Rc from the collapse; we thus expect 
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rA ~ Rc ~ ^0 and hence /? « 1. Notice that in the hmit r ^ 0, a ^ 0, and the function takes 
on a hmiting form h*{x; a 0) = a/tt/S e^ Erf[-y/x] — \/2x(l + x). 

Figure 2 shows the density profiles for a representative cluster with an anisotropic velocity 
distribution, where the parameters are taken to be ipQ = 10, A = 1/2, and (3 = 1. The dashed 
curve shows the stellar component, the dotted curve shows the gaseous component, and the solid 
curve shows the total density; the dot-dashed curve shows the reduced stellar density profile after 
gas removal (see the following section). For this class of models, the velocity anisotropy leads to 
a more extended stellar component compared to the isotropic models of the previous section. For 
the same values of -00 and A, the stellar component extends farther into the gaseous halo and the 
mass fraction of stars is correspondingly smaller. For this typical model using ipQ = 10 and A = 
1/2, for example, the stellar fraction is 27% (compared with 44% for the isotropic case /3 = 0). The 
variation of the density profiles with varying depths of the gravitational potential {il^o) and varying 
amounts of gas (A) are similar to those shown in Figure 1 for isotropic models. 



For the next stage of this calculation, we let all of the gas be removed from the cluster on a 
short time scale. We thus assume that the distribution function of the stars does not have time to 
adjust and hence all of the stars suddenly find themselves in a new environment with a less deep 
gravitational potential well, but they initially retain the same distribution function, as given by 
equation [4] or [12] evaluated using the old values of the potential. In nature, the removal of gas 
does not take place instantaneously, however, so this calculation represents a limiting case of the 
physical problem. 

The energy required to remove gas from a cluster can be easily obtained. The binding energy 
of a cluster is given hy E = kGM'^/R, where A; is a constant of order unity that depends on the 
cluster shape. In terms of representative values, we can write 



Since a supernova explosion has a typical energy rating of Esn ~ 10 erg, one supernova provides 
enough energy to have a devastating impact on the gaseous component of a cluster. However, 
radiative cooling is very efficient in molecular clouds and can dispose of much of the injected energy 
on short time scales (Wheeler et al. 1980; Franco et al. 1994). In addition, gas dispersal takes place 
efficiently even in the absence of supernovae (Palla & Stabler 2000). As another mechanism, winds 
from young stellar objects inject a large amount of mechanical luminosity into the surrounding 
medium, with a typical scaling law of ^mech ~ O.OIL* (e.g., Lada 1985). According to this relation, 
a young stellar object with = 3 imparts an energy AE ~ 10^^ erg during a time interval 
~0.3 Myr. The outflows from YSOs thus impart enough energy to facilitate gas removal. At the 
high densities of molecular clouds, however, gas removal by photoionization and photodissociation 
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REMOVAL OF THE GASEOUS COMPONENT 




(23) 
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(AV'g)rT = A(l - — ata<T) ■ (25) 



can be a more effective mecfianism (Diaz-Miller et al. 1998). In any case, tfie energy to remove gas 

from clusters is easily obtained, but the requisite momentum transfer is more difficult to achieve. 
Although the actual gas dispersal mechanism remains unclear, gas is observed to dissipate in less 
than 10'' years, even for aggregates containing no massive stars (Palla & Stahler 2000). 

In the models of this paper, we calculate the fraction of stars remaining in a cluster (after 
gas removal) as a function of star formation efficiency. We note, however, that the star formation 
efficiency is coupled to the gas removal process. If gas is removed slowly, then the transformation of 
gas into stars operates over a longer period of time and the star formation efficiency can be larger 
(although the opposite is not necessarily true). 

Against this background, we now proceed with the calculation. The immediate change A^g 
in the gravitational potential - that just due to gas removal - is given by 

AV'g(e) = (AV'g),^ + A [^^atan^T - ^atan^ + ^ ln(^^^^)] , (24) 

where the change in potential (Ax/jg)rj, at the outer cluster boundary can also be evaluated to 
obtain ^ 

^' 

This form for the outer boundary condition assumes that the gaseous halo is embedded within a 
much larger gas cloud with the properties of an isothermal sphere. This assumption thus specifies 
the constant $o = 0- We can combine the above two expressions and simplify the result to obtain 

A^.iO = A [l - latan^ + i ln(^^^)] ■ (26) 

The physical change in the potential is thus given by A\I'g = Aipg{^). 

With the gas removed from the cluster, the stars find themselves in a new (shallower) potential 
and hence some stars must leave the system. The gravitational potential changes because of both 
gas removal and the removal of the high velocity stars. At a given location within the cluster, the 
original distribution of stars had a range of possible velocities given by 

0<v'^ <2^ . (27) 

In the new configuration, without the gaseous component, the stars that remain bound to the 
cluster must have velocities in the range 

< < 2(* - A^') = v^"^ , (28) 

where A'l' includes the change in potential due to stars leaving the system, and where we have 
defined a maximum velocity in the final equality. We thus need to solve for the new potential, 
which we can write as ^ = v^^ /2a'^ . 
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3.1. Reduced density profiles for Isotropic models 

For our model with isotropic velocity dispersion, the density of stars remaining is given in 
terms of (p (or, equivalently, Vm.) according to the relation 

/9*(bound)= / 4Trv'^dvf{e) = pi e'^Erf(^) - 2V^(e'^"'^ + -0) . (29) 

In the above expression, the old potential ip no longer has the physical meaning of the gravitational 
potential. Instead, V'(0 is simply a known function which determines the distribution function of 
stars, only a fraction of which will remain in the cluster. The fraction that remains depends on the 
new potential ^(^) which must be determined by solving the following new version of the Poisson 
equation, 



r'i(e^) = --^U''mV^)-2VH^{e^-^ + 'U)] , (30) 



2 



3 



with i^iC) known function. To find the new potential, and hence the new density, we must also 
specify the boundary condition at the cluster center ^ = 0. Here we invoke the condition 

m = m - Av^g(o) = m - ^ hi + ^t^) , (31) 

which accounts for the change in potential due to gas loss. 

With the problem now completely specified, we can solve for the original equilibrium structure 
of the cluster as a function of the potential well depth and the gas parameter A. We can then 
integrate equation [30] to find the density profile for the stars remaining after the gas is removed. 

The mass fractions are listed in Table 1 for the representative case of ipQ = 10. Notice how the 
results differ from those predicted by equation [1]. For A = 1/2, for example, the total gas mass is 
somewhat greater than the original mass in stars. After the gas is removed, nearly 73% of the stellar 
component remains, whereas the virial argument predicts that none of the stars remain bound. 

The difference between these results lies in two key features of this present treatment: [A] By 
formulating the problem in terms of the distribution function /(e) for the stars, we take into account 
the low velocity tail of the distribution (which always tends to remain bound). [B] The gaseous 
halo in our models is tied to the background molecular cloud, and not directly tied to the stellar 
component of the cluster (although the gravitational interaction is included). This complication 
allows the gas density distribution to have a markedly different form from that of the stars. In 
particular, the stellar component exhibits a relatively steep fall off near the edge of the cluster, 
whereas the gas density continues to follow its usual power-law distribution (see Figure 1). As a 
result, the mass of the gaseous component is concentrated more towards the outer portion of the 
cluster (compared to the stellar component) and gas removal has a less destructive effect on stellar 
population of the cluster. 

Figure 3 shows the fraction JT^ of stars remaining after gas removal as a function of the star 
formation efficiency, defined as e* = / p^dV/M. Notice that the shape of these curves is markedly 
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different from the singular behavior predicted by equation [1]. As a very rough approximation, 
the entire collection of curves can be described by the simple function = 2e* — e*, as shown 
by the dashed curve. Much more accurate fits to individual models can be obtained. For the 
t/jQ = 10 model, for example, the dotted curve shows a fit using a function = (2?— e^)^/^, where 
?= (lOe* — l)/9 is a stretched variable. 

The clusters produced by this rapid gas removal have a truncated distribution function (by 
construction). Prom this state, the cluster tends to evolve toward a new configuration with a 
smooth distribution function. Because this subsequent evolution takes place through dynamical 
interactions of the constituent stars, the time scale for readjustment is longer than the gas removal 
time (but is nonetheless shorter than the total dynamical evolution time scale). This longer term 
evolution poses an interesting problem for future work. 



3.2. Reduced density profiles for Anisotropic models 

For the anisotropic models developed in the previous section, the new density profile, after the 
gas is removed from the cluster, can be written in the form 

P* = Piy||e^i^p(0,a) + ^atan(aV2^) [e^-^ + + -L] | . (32) 

Keep in mind that V' is the original potential, before gas removal, and is the physical potential 
as determined by the solution to the Poisson equation. 

As in the previous case, we can find the equilibrium structure of the cluster in the presence 
of a gaseous component, and then use the above fornuilation to find the density profile of the 
stars remaining after gas is removed. Compared with the isotropic models of the previous section, 
the anisotropy of the velocity distribution allows a greater fraction of the stars to remain after 
gas removal. For a collection of models with varying i/jq, Figure 4 shows the fraction of stars 
remaining after gas removal as a function of the star formation efficiency e* (the solid curves) . The 
dashed curve shows the function J^^ = (2e* — e^)^/^, which provides a rough fit to the family of 
curves as shown. As illustrated here, the anisotropic models preserve an even greater fraction of 
their stars than the isotropic models considered previously (compare Figure 3 and Figure 4). 



4. COMPARISON WITH OBSERVATIONS 

The Trapezium cluster is a nearby young cluster in Orion and is often used an example of 
a forming bound cluster (e.g., Hillenbrand & Hartmann 1998; McCaughrean & Stauffer 1994; 
Hillenbrand 1997). Since this cluster is relatively young (the stars have a quoted mean age of ~ 0.8 
Myr), it provides a good point of comparison for the models of newly formed clusters. The total 
mass within the central 2.06 pc region is estimated to be in the range 4500 - 4800 M© (Hillenbrand 
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fe Hartmann 1998); these same authors estimate that the core radius ro ~ 0.2 pc and find a total 

number of stars A^* ~ 2300 within the central r = 2.06 pc region. Within this same volume, the 
total mass of the cluster is larger than that of the stellar component by a factor of about two, i.e., 
stars constitute 40-50% of the cluster mass within the inner 2 pc. 

The characteristics of this cluster are roughly consistent with the cluster formation scenario 
outlined in the Appendix. If we use the isothermal version of the model, for example, we would 
need an effective sound speed of a 2.3 km/s to form a 4800 M© cluster in 1.7 Myr; these same 
values imply Too ~ 2 pc. In order to obtain a centrifugal radius (and a core radius) of Rc = tq = 0.2 
pc, the initial rotation rate of the core would have to be Hi 0.4, which is a reasonable value for 
molecular clouds. Similarly, for the logatropic model, a pressure scale of Pq ~ 6.6 xlO~^'^ dyn 
cm~^ would imply r^o = 2 pc and a formation time of 2.9 Myr. To obtain the same core radius 
of ro = Rc = 0.2 pc, the required rotation rate is again f^i fs 0.4. We thus argue that inside-out 
collapse models can be made consistent with this particular observed young cluster. 

Hillenbrand & Hartmann (1998) have already used lowered isothermal models to fit to the 
stellar component of the Trapezium cluster data and find good agreement using a central potential 
V'o = 9. This model does not include the gaseous component or the possibility of anisotropy in the 
velocity distribution; however, including these additional parameters can only make the fit better. 
We can obtain a reasonable model fit to this same cluster using -00 ~ 10. For an isotropic model 
with A = 1/2, the total mass in mass is comparable to the mass in stars. If the gas is removed 
from such a cluster, nearly 73% of the stellar mass remains bound, with the remainder escaping to 
fill the field with stars. These models extend out to a truncation radius tt which is about 100 core 
radii and thus extends many times farther out than the observed portion of the cluster. In order 
to have the gas mass and the stellar mass nearly equal in the (observed) central 2 pc region, we 
need a larger gas parameter of A ~ 1.5; for this case, only about 48% of the stars remain bound 
after the gas leaves the system. After gas removal, this simple model thus predicts that the central 
region of the Trapezium cluster will remain a bound entity with approximately ~ 1000 stars. 
If the velocity distribution is not isotropic, but instead becomes more radial in the outer cluster, 
then an even larger fraction of stars could be retained after gas removal. 

In the existing literature, there are few predictions regarding the fate of the Trapezium for 
comparison. An older study using only 16 stars indicates a half-life of only 10^ years (Allen & 
Poveda 1975). On the other hand, a recent numerical study (Kroupa 2000) obtains results in good 
agreement with the semi-analytic models of this paper. 

5. SUMMARY and DISCUSSION 

In this paper, we have used results from star formation theory and stellar dynamics to develop a 

working model for the early evolution of open star clusters. More specifically, we have constructed a 
sequence of equilibrium cluster models, including both isotropic and anisotropic velocity dispersions. 



- 12 - 



These models represent a straightforward generahzation of previous work to include a gaseous 
component (see Figures 1 and 2). We have then considered the effects of gas removal from these 
clusters. 

When a cluster loses its gaseous component, the subsequent evolution of the stellar component 
is largely determined by the distribution function /(e,L). The cluster always contains some stars 
with low velocities (low kinetic energy), and hence at least some stars always remain behind. 
Because the distribution function depends sensitively on the velocity, however, the number of stars 
left behind is sensitive to the shape of the distribution. In this class of models, relatively more stars 
remain bound to the cluster than suggested by virial arguments. For example, the model most 
applicable to young open star clusters, with anisotropy parameter (3 = 1 and with equal mass in 
gas and stars, initially retains 93% of its stars after the gaseous component is removed. In general, 
the fraction ^*(e*) of stars left behind after gas removal is a smoothly varying function of star 
formation efficiency (see Figures 3 and 4). As a crude approximation spanning the entire space 
of solutions for isotropic models, the fraction can be characterized by a elementary function of 
the form J^^ = 2e* — e^. More accurate fits can be obtained for specific cluster models. Similarly, 
anisotropic models with P = 1 can be described by the function JT^ = (2e* — e^)^/^. 

Along the way to the above results, we have developed a simple model of cluster formation 
(described in the Appendix). If clusters form out of the collapse of a large molecular cloud structure, 
the collapse flow can be described by scaled-up versions of the infall collapse solutions for individual 
star formation events. The time scale for individual stars to form is much shorter than the time 
required for the cluster to form, and the dynamics on the larger size scale of the cluster become 
purely ballistic. The central portion of the forming cluster, at radii smaller than the centrifugal 
radius Rc ~ 0.1 pc, has time to dynamically relax and tends to exhibit an isotropic velocity 
dispersion; the centrifugal radius also sets the core radius of the cluster. In the outer portion of 
the forming cluster, r > Rc, the velocity dispersion becomes anisotropic and nearly radial. 

This paper represents a modest step toward a unified picture of cluster formation, early evolu- 
tion, and dispersal. We have included a gaseous component in the construction of semi-analytical 
equilibrium models of clusters and have considered the initial effects of gas removal. The subse- 
quent development of the resulting stellar systems should be studied next. More detailed models of 
cluster formation constitute another fruitful area for further work. With a more definitive picture 
of cluster formation in hand, the distribution function of the stars in the earliest phases can be more 
precisely specified and the fraction of stars that remain after gas removal can be more accurately 
determined. Another important complication is to include a different distribution function for stars 
of different masses. Once we understand these issues of cluster formation and early evolution, we 
can then determine their effects on the star formation process. 

We would like to thank Gus Evrard, Charlie Lada, Greg Laughlin, Phil Myers, Doug Richstone, 
and Frank Shu for useful discussions. This work was supported by funding from The University of 
Michigan and by bridging support from NASA Grant No. 5-2869. 
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A. A Cluster Formation Model 

In this Appendix, we present a simple model for cluster formation. For robust bound clusters 
forming within molecular clouds, we expect the original proto-cluster structure to collapse as a 
whole. To obtain a mathematical description of this collapse, we consider the flow that eventually 
produces a cluster to be a scaled up version of the collapse flows that have been studied previously 
for single star formation (e.g., Shu 1977; Terebey, Shu, &: Cassen 1984; Jijina &; Adams 1996). This 
approach should capture the basic essence of the collapse problem. In this case, the collapse of 
a molecular cloud region proceeds from inside-out. The central portion of the flow approaches a 
ballistic (pressure-free) form and helps determine the velocity distribution of the forming cluster. 
Even for infalling gas, the inner limit of the collapse flow always approaches pressure-free conditions; 
infalling stars are manifestly ballistic. The time scale for individual star formation events is ~ 
10^ years (Myers &: Fuller 1993; Adams & Fatuzzo 1996), whereas the time scale for the entire 
cluster to form is somewhat longer, 1-3 Myr. We thus expect most of the stars to form while the 
overall collapse of the cluster is still taking place. In apparent support of this picture, observations 
suggest that cluster formation takes place within only about one sound crossing time of the system 
(Elmegreen 2000). 

For a given gravitational potential, we find the orbital solutions for stars (or gas parcels) falling 
towards the cluster center. In the standard infall calculation, the inner solution is derived using 
the gravitational potential of a point source. Since this potential is spherically symmetric, angular 
momentum is conserved and the motion is confined to a plane described by the coordinates (r, cp); 
the radius r is the same in both the plane and the original spherical coordinates. The angular 
coordinate (f) in the plane is related to the angle in spherical coordinates by the relation cos0 = 
cos 9 I cosOq, where is the angle of the asymptotically radial streamline (see below). For zero 
energy orbits, the equations of motion imply a cubic orbit solution, 

l-^ = C(l-/.g). (Al) 

Here, the trajectory that is currently passing through the position given by ^ and jU = cos 6 initially 
made the angle (a*o = cos^o) with respect to the rotation axis. The quantity ^ is defined by 

where joo is the specific angular momentum of parcels of gas currently arriving at the cluster center 
along the equatorial plane. The second equality defines a centrifugal radius Rc- We assume that 
the initial cloud is uniformly rotating at a constant rotation rate O, so that joo = f^r^, where Too 
is the starting radius of the material that is arriving at the center at a given time. 

To evaluate the radii Rc and Too, we invert the mass distribution of the initial state. For an 
isothermal cloud, the mass profile and the centrifugal radius are given by 

^ 2aV GM , „ n^G^M^ 

M{r) = -^, r^ = ^, and Rc = , (A3) 
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where a is the isothermal sound speed. For this case, the infall collapse solution (Shu 1977) indicates 
that the flow exhibits a well defined mass infall rate M = rriQa? /G, where mg ~ 0.975. To form a 
cluster with mass M = IOOOMq in 1 Myr, for example, the required effective sound speed is a ~ 
1.63 km/s. With this central mass, the region that originally filled a volume of radius Too ~ 0.81 
pc has fallen to the center (within Rc)- The size of the collapsing region is about twice as large, 
th = at 1.66 pc, and contains a total mass M 2000 Mq. The centrifugal radius Rc ~ 0.1 pc 
for f2 = 1 km s"^ pc~^ (a typical observed value; Goodman et al. 1993). This centrifugal radius 
is comparable to the expected core radius rcore of a newly formed cluster, and we make the rough 
identification rcore ^ Rc- 

We also consider initial states for non-isothermal conditions. Molecular linewidths often show 
a substantial nonthermal component with a density dependence of the form Av oc p~^^'^ (e.g., 
Larson 1985; Myers & Fuller 1992). If we use a "logatropic" equation of state P = Pologp/po to 
describe such a fluid (Lizano Sz Shu 1989; Jijina & Adams 1996; McLaughlin & Pudritz 1997; Galli 
et al. 1999), the equilibrium mass distribution and the corresponding radii r^o and Rc are given 

by 



M(r) = Too = — , and Rc = , (A4) 



G 



27rPo 



G 



-1/4 , _ n'^M 



27rPn 



where Pq is the pressure scale that determines the amount of nonthermal support in the cloud. In 
this case, the mass infall rate is time dependent. The total mass M{t) that falls to the center of the 
flow during a time t is given by M = mot^{2TrGPo)^^^ /16G, where mo ~ 0.0302. During logatropic 
collapse, most of the mass in the original cloud region is still on the way down, rather than at the 
center. If the cluster encompasses the entire collapsing region, the total mass is about 30 times 
that of the central core. A typical scenario would thus have 100 Mq in the central core and 3000 
Mq still falling inwards. For this case, the required pressure scale Pq ~ 8.9 x 10~^'^ dyne/cm^. To 
obtain a centrifugal radius Rc ~ 0.1 pc ~ rcore as before, we need $7 3 km pc~^. 

Both the isothermal and the logatropic models can produce a cluster in ~ 1 Myr using rea- 
sonable values for the input parameters {CI and a or Pq). In both cases, the effective transport 
speed required to support the initial cloud is comparable to the velocity dispersion of the resulting 
cluster. For the range of parameters discussed above, this velocity scale is 1-2 km/s and is roughly 
consistent with the velocity dispersion of observed open clusters. 

For this chistcr formation scenario to be consistent with the current paradigm of star formation 
(e.g., Shu, Adams, &: Lizano 1987), the cluster environment cannot greatly disrupt the collapse of 
smaller regions that produce individual stars. To fix ideas, wc assume that the large scale collapse 
of the cluster is determined by logatropic conditions (eq. [A4] ) and the collapse of individual star 
forming regions is given by isothermal conditions (eq. [A3]). To produce a 1 M© star, for example, 
the radial extent of the initial pre-collapse region is roo = GM/la^ 0.02 pc. In the central region 
(r < 1 pc) that eventually becomes the cluster, the individual star forming sites do not greatly 
interfere with each other as long as the number of stars docs not exceed N ~ (1/0.02)'^ 10^. 
Similarly, the collapse of an individual star proceeds largely independent of the tidal forces. For a 
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cluster of mass M^ust ~ 1000 Mq and size -Rciust ~ 1 PC, the tidal radius ~ 0.1 pc, which is 
much larger than the size of an individual infall region (~0.02 pc). 

We expect that the stars in the newly formed cluster retain some dynamical memory of the 
velocity distribution of the collapse flow. In this flow, streamlines entering the central region do 
not cross each other. As long as the infall time is longer than the time scale for individual star 
formation events, the core regions that produce stars will not have a chance to interact. [As an 
aside, note that this cluster formation scenario has an initial transient phase (the first 10^ yr) in 
which the infall of the cluster takes place faster than individual stars can form. Protostars entering 
the central region during this initial time period thus have an opportunity to interact and merge; 
this activity may contribute to the production of more massive stars in the cluster center.] 

Given the orbital solution, we can find the velocity fields for the collapse flow: 

.r^-[^y\2-ai-,l)y^\ (A5) 



-o={^] ' \'^A^^o-^^')C} , (A6) 




1/2 r _ o ^ 1/2 

2 .,2\ 



ii-f^Dii-f^'r'/'c'^'. (A7) 



Since /i, and fiQ are related through the orbit equation [Al], the velocity field is completely 
determined for any position {r,6). With the velocity field specified, we can find the anisotropy 
in the flow as a function of radius. We define the angular average for both the perpendicular 
component of the velocity field and the radial component 

(^i) = {ve +vl) = ^CIv and {vD = ^ [2 - C4] , (A8) 



where the integral 1^ is given by 



Jo 



/ d/x(l-//^). (A9) 



To evaluate the integral ly, we change the integration variable from /i to fj,o and change the 
lower end of the range of integration from to a critical value /xc. This difference arises because 
streamlines from all initial angles cannot fall to arbitrarily small radii. For large radii, streamlines 
from all values of hq are represented. Inside the centrifugal barrier Rc, however, only streamlines 
originating preferentially from the poles reach these smaller radii. The last streamline that reaches 
a given radius is determined by iJ,c- Evaluating the integral we find 

4 = (|-^C)(i-Mc) + M4 + ^i). (AlO) 
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This expression simplifies in the inner and outer regimes. For large radii, r Rc, fJ-c 0; ^^nd 
/„ (2/3)(l-2C/5). In the opposite limit of small radn, r < Rc, He ^ 1-1/2C, and /„ 8/15C. 

In the context of cluster formation, we evaluate the anisotropy of the flow in the outer regime 
and assume that the cluster retains some memory of its initial velocity distribution. Combining 
equations [A8 - AlO], we find 

^ = _ _ c (1 - 2C/5) 

To leading order, we thus obtain Tly ^ C, = Rc/r. For large radii the velocities become nearly 
radial, whereas for small radii the velocities become more isotropic. Scattering of the newly formed 
stars will be relatively efficient at small radii and can drive the velocity dispersion even further 
towards isotropy. More specifically, the relaxation time can be written in the form trclax ~ 14 Myr 
(r/lpc)^ (roo/lpc)~^, where we have assumed a typical robust cluster with N = 1000 stars. For 
small radii r < 0.1 pc Rc-, the relaxation time is less than about 0.2 Myr. As a result, the region 
within the centrifugal radius will experience several relaxation times during the expected time scale 
for the cluster to form. 

In this picture, clusters form with a nearly isotropic velocity dispersion on the inside and 
a highly radial velocity dispersion on the outside, with the centrifugal radius Rc providing the 
boundary between the two regimes. Furthermore, the core radius of the cluster is determined by 
ro ~ Rc ~ 0.1 pc (for typical initial conditions). 
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Table 1. Parameters for Isotropic models with i/;o = 10 



A 












0.00 


224 


126 


0.00 


126 


1.000 


0.25 


150 


84.9 


51.7 


73.4 


0.865 


0.50 


106 


58.9 


73.0 


42.8 


0.727 


0.75 


78.1 


42.7 


80.1 


26.9 


0.631 


1.00 


59.6 


32.5 


81.0 


18.5 


0.569 


1.25 


47.1 


25.8 


79.4 


13.6 


0.526 


1.50 


38.3 


21.2 


77.0 


10.4 


0.489 


1.75 


32.0 


17.9 


74.3 


8.09 


0.451 


2.00 


27.2 


15.5 


71.8 


6.36 


0.412 


2.25 


23.6 


13.6 


69.4 


4.99 


0.369 


2.50 


20.8 


12.0 


67.2 


3.88 


0.323 


2.75 


18.5 


10.8 


65.3 


2.96 


0.274 


3.00 


16.7 


9.76 


63.5 


2.19 


0.225 


3.25 


15.1 


8.90 


61.9 


1.55 


0.174 


3.50 


13.9 


8.16 


60.4 


1.00 


0.123 


3.75 


12.8 


7.53 


59.1 


0.52 


0.069 


4.00 


11.9 


6.97 


57.9 


0.09 


0.013 
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FIGURE CAPTIONS 

Fig. 1. — Equilibrium density profiles for clusters with isotropic velocity distributions. The density 
is scaled to the central density po of the stellar component and the radius is given in units of the 
corresponding King radius, i.e., ^ = r/ro. (a) A typical density profile with i/^q = 10 and A = 1/2 

(44% stars and 56% gas). The dashed curve shows the stellar component, the dotted curve shows 
the gaseous component, and the solid curve shows the total density. The dot-dashed curve shows 
the stellar component remaining after gas is removed from the system, (b) Collection of density 
profiles with varying depths of the gravitational potential and fixed gaseous halo with A = 1/2. 
The various curves use V'o = 6, 8, 10, 12, and 14. As the depth V'o increases, the outer radius of 
the stellar component of the cluster grows accordingly (see also Tabic 1). (c) Collection of density 
profiles for ipQ = 10 and varying amounts of gas, determined by A = 0, 0.5, 1.0, 1.5, and 2.0. As A 
increases, the outer radius of the cluster decreases. 

Fig. 2. — Equilibrium density profile for a cluster with an anisotropic velocity distribution. The 
density and radius have the same units as in Figure 1. This model uses ipQ = 10, A = 1/2, and (3=1. 
The dashed curve shows the stellar component, the dotted curve shows the gaseous component, 
and the solid curve shows the total. The dot-dashed curve shows the stellar component remaining 
after gas is removed from the cluster. 

Fig. 3. — Fraction .F* of stars remaining bound to a cluster after gas removal for systems with 
isotropic velocity distributions. The solid curves show the fraction JR^ as a function of star formation 
efficiency e* for varying ^/jq (here, tpQ = 6, 8, 10, 12, and 14). The dashed curve shows the function 
= 2 e* — e^, which provides a rough approximation to the entire family of curves. The dotted 
curve shows a more accurate fit (to the ?/'o = 10 model) using the function = (2?— e^)^/^, where 
?= (lOe* — l)/9 is a stretched variable. 

Fig. 4. — Fraction J^^ of stars remaining bound to a cluster after gas removal for systems with 
anisotropic velocity distributions. Curves show the fraction as a function of star formation 
efficiency e* for varying ipQ (here, Vo = 8, 10, and 12). This family of curves was calculated using 
velocity anisotropy (5=1. The dashed curve shows the function J^^ = (2e* — e^)^/^, which provides 
a rough fit to the family of curves as shown. 
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